TP Clustering - Correction

Partie 5 : Clustering des données Vélib

library(funFEM)
library(reshape2)
library(ggplot2)
library(gridExtra)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(mclust)
library(stringr)
library(cluster)
library(clusterSim)
library(ggmap)

1 Etude préliminaire

1.1 Préparation des données

Cet ensemble de données contient des données provenant du système de partage de vélos de Paris, appelé Vélib. Les données sont des profils de chargement des stations de vélos sur une semaine. Les données ont été collectées toutes les heures pendant la période du dimanche 1er septembre au dimanche 7 septembre 2014.

On charge ici les données Velib disponibles dans la librairie funFEM.

library(funFEM)
data(velib)

Ces données contiennent

• data : les profils de chargement (nb de vélos disponibles / nb de points d’attache) des 1189 stations à 181 points de temps.

• position : la longitude et la latitude des 1189 stations de vélos.

• dates : les dates de téléchargement.

• bonus : indique si la station est sur une colline (bonus = 1).

• names : les noms des stations.

Afin d’avoir des semaines complètes, nous allons supprimer les 13 premières colonnes.

Velibdata <- velib$data[, -c(1:13)]
colnames(Velibdata) <- velib$dates[-c(1:13)]

On controle le nombre de jours et d’heures dans le jeu de données

day <- str_sub(colnames(Velibdata), 1, 3)
day <- factor(day, levels = c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim"))
table(day)

day
Lun Mar Mer Jeu Ven Sam Dim 
 24  24  24  24  24  24  24 

hour <- str_sub(colnames(Velibdata), 5, 6)
table(hour)

hour
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
 7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 

timeTick <- 24 * (0:6)

On stocke les informations sur les stations dans le dataFrame suivant :

station <- data.frame(nom = velib$names, velib$position, colline = velib$bonus)

On ajoute des informations complémentaires sur les stations :

# Les gares
gare <- rep(0, nrow(station))
gare[which(str_detect(station$nom, "GARE") == T)] <- 1

# Lieux touristiques (autour de Tour eiffel, Panthéon, Arc de Triomphe, Champs
# Elysées, Louvre, Notre Dame, Invalides, Trocadéro, Sacré Coeur, Montmartre)
tourisme <- rep(0, nrow(station))
TI <- c(209, 504, 897, 102, 518, 357, 256, 301, 992, 627, 579, 241, 646, 13, 283,
    941, 976, 1189, 22, 633, 853, 973, 762, 195, 71, 131, 912, 302, 303, 1068, 774,
    384, 307, 326, 436, 789, 250, 755, 461, 439, 918, 1065, 99, 151, 313, 871)
tourisme[TI] <- 1

# Les mairies
mairie <- rep(0, nrow(station))
mairie[which(str_detect(station$nom, "MAIRIE") == T)] <- 1
mairie[which(str_detect(station$nom, "HOTEL") == T)] <- 1

# Bois / Parc / Square
parc <- rep(0, nrow(station))
parc[which(str_detect(station$nom, "BOIS DE") == T)] <- 1
parc[which(str_detect(station$nom, "PARC") == T)] <- 1
parc[which(str_detect(station$nom, "SQUARE") == T)] <- 1
parc[which(str_detect(station$nom, "CANAL") == T)] <- 1

# Porte
porte <- rep(0, nrow(station))
porte[which(str_detect(station$nom, "PORTE") == T)] <- 1

station <- data.frame(station, gare = gare, tourisme = tourisme, mairie = mairie,
    parc = parc, porte = porte)

1.2 Quelques fonctions auxiliaires pour exploiter les résultats

Fonction pour tracer les profils (charge de chaque station / loading) en fonction des jours et heures pour les stations choisies.

plotprofils <- function(Velibdata, station, numstation, plotsem = TRUE) {
    # Velibdata = les données de velib initiales station = ens. des
    # informations sur les stations numstation = vecteur contenant le numéro de
    # ligne des stations dont on souhaite tracer les profils

    library(reshape2)
    library(ggplot2)
    Dataaux <- data.frame(id.s = station$nom, Velibdata)
    Aux <- melt(Dataaux[numstation, ], id = c("id.s"))
    Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
        5, 6)))
    Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
        "Ven"))

    if (plotsem == TRUE) {
        g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
            ncol = 2) + xlab("Time") + ylab("Loading") + scale_x_continuous(breaks = seq(0,
            24, 4)) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
            labs(col = "Classif.") + scale_x_continuous(breaks = c(0, 7, 9, 12, 14,
            20, 24))
    } else {
        v <- rep(0, nrow(Aux))
        for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
            levels(Aux$day)[j])] + (24 * (j - 1))
        Aux <- data.frame(Aux, v = v)
        rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
            ymax = rep(1, 7), color = rainbow(n = 7))
        g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
            geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
                fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
            ylab("Loading") + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
            labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
            7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
    }
    return(g)
}

Exemple d’utilisation :

plotprofils(Velibdata, station, numstation = c(813, 655, 468))

plotprofils(Velibdata, station, numstation = c(813, 655, 468), plotsem = FALSE)

Se donnant une classification, fonction pour tracer les profils moyens par classe :

plotprofilsmoy <- function(Velibdata, clustering, plotsem = TRUE) {
    library(reshape2)
    library(ggplot2)

    Dataaux <- matrix(0, nrow = max(clustering), ncol = ncol(Velibdata))
    for (k in 1:max(clustering)) {
        Dataaux[k, ] <- apply(Velibdata[which(clustering == k), ], 2, mean)
    }
    colnames(Dataaux) <- colnames(Velibdata)
    Dataaux <- data.frame(id.s = as.factor(1:max(clustering)), Dataaux)
    Aux <- melt(Dataaux, id = c("id.s"))
    Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
        5, 6)))
    Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
        "Ven"))

    if (plotsem == TRUE) {
        g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
            ncol = 2) + xlab("Time") + ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") +
            theme(axis.title.x = element_blank()) + labs(col = "Classif.") + scale_x_continuous(breaks = c(0,
            7, 9, 12, 14, 20, 24))
    } else {
        v <- rep(0, nrow(Aux))
        for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
            levels(Aux$day)[j])] + (24 * (j - 1))
        Aux <- data.frame(Aux, v = v)
        rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
            ymax = rep(1, 7), color = rainbow(n = 7))
        g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
            geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
                fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
            ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
            labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
            7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
    }
    return(g)
}

Exemple d’utilisation avec les stations sur colline par rapport aux autres :

plotprofilsmoy(Velibdata, clustering = station$colline + 1, plotsem = TRUE)

plotprofilsmoy(Velibdata, clustering = station$colline + 1, plotsem = FALSE)

Fonction auxiliaire pour comparer une classification avec les informations auxiliaires disponibles pour les stations

stationcaract <- function(station, clustering) {
    library(dplyr)
    # Dataaux<-data.frame(id.s=c(1:nrow(Data)),station) aux <- cbind(Dataaux,
    # clust)
    aux <- cbind(station, clustering)
    aux.long <- melt(data.frame(lapply(aux[, -c(2:3)], as.character)), stringsAsFactors = FALSE,
        id = c("nom", "clustering"), factorsAsStrings = T)
    # Effectifs
    aux.long.q <- aux.long %>%
        group_by(clustering, variable, value) %>%
        mutate(count = n_distinct(nom)) %>%
        distinct(clustering, variable, value, count)
    # avec fréquences
    aux.long.p <- aux.long.q %>%
        group_by(clustering, variable) %>%
        mutate(perc = count/sum(count)) %>%
        arrange(clustering)

    gaux <- ggplot(data = aux.long.p, aes(x = clustering, y = perc, fill = value)) +
        geom_bar(stat = "identity") + facet_wrap(~variable)

    return(gaux)
}

1.3 Quelques statistiques descriptives

On reprend ici rapidement quelques éléments d’analyse descriptive vus précédemment dans l’UF.

• Tracé des corrélations des heures par jour, toutes les stations confondues :

library(corrplot)
timeTickAux <- c(timeTick, 168)
for (k in 1:7) corrplot(cor(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k + 1])]),
    method = "ellipse")

• Comportement moyen par jour :

velibmeanday <- matrix(0, nrow = nrow(Velibdata), ncol = 7)
for (k in 1:7) {
    velibmeanday[, k] <- apply(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k +
        1])], 1, mean)
}
colnames(velibmeanday) <- c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim")
ggplot(melt(velibmeanday), aes(x = Var2, y = value)) + geom_boxplot() + xlab("") +
    ylab("Loading")

• Comportement moyen par heure et par jour :

velibHour <- data.frame(value = colMeans(Velibdata), day = day, hour = as.numeric(hour))
ggplot(velibHour, aes(x = hour, y = value, col = day)) + geom_line() + geom_point() +
    ylab("Loading") + ggtitle("Hourly loading, averaged over all stations")

• Quelques visualisations des stations sur le plan 2D de Paris :

library(ggmap)
# Average loading for each station
load <- rowMeans(Velibdata)
qmplot(longitude, latitude, data = station, color = load) + scale_colour_gradient(high = "red",
    low = "yellow")

# scale_colour_gradient(high='blue',low='white')

# Loading moyen à l'heure t (le premier temps est 0h)
t <- 0
loadheure <- rowMeans(Velibdata[, seq(from = (t + 1), by = 24, length = 7)])

qmplot(longitude, latitude, data = station, color = loadheure) + scale_colour_gradient(high = "red",
    low = "yellow")

• Visualisation des stations selon les informations auxiliaires sur les stations :

qmplot(longitude, latitude, data = station, color = as.factor(station$colline))

qmplot(longitude, latitude, data = station, color = as.factor(station$tourisme))

1.4 Analyse en composantes principales

Question : Faites une ACP des stations et analysez les résultats.

respca <- PCA(.......)
fviz_eig(....)
fviz_pca_var(......)
fviz_pca_ind(........)

Correction :

library(FactoMineR)
library(factoextra)
respca <- PCA(Velibdata, scale.unit = T, graph = F, ncp = 15)

# Variance expliquée
fviz_eig(respca)

ggplot(data.frame(comp = 1:20, cumul = respca$eig[1:20, 3]), aes(x = comp, y = cumul)) +
    geom_line() + geom_point()

# Visualisation des variables
fviz_pca_var(respca, axes = c(1, 2), geom = c("arrow"))

fviz_pca_var(respca, axes = c(2, 3), geom = c("arrow"))

for (k in 1:7) corrplot(respca$var$cor[(timeTickAux[k] + 1):(timeTickAux[k + 1]),
    1:3], method = "ellipse")

# Projection des individus
fviz_pca_ind(respca, axes = c(1, 2), geom = c("point"), habillage = as.factor(station$colline))

fviz_pca_ind(respca, axes = c(1, 3), geom = c("point"))

nbdim <- 4
aux <- melt(respca$var$coord[, 1:nbdim])
aux$Var1 <- str_sub(aux$Var1, 1, 3)
aux <- data.frame(aux, time = rep(c(0:167), nbdim))

rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(min(aux$value),
    7), ymax = rep(max(aux$value), 7), color = rainbow(n = 7))

ggplot(data = aux) + geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin,
    ymax = ymax, fill = color), alpha = 0.3) + geom_line(data = aux, aes(x = time,
    y = value)) + facet_wrap(~Var2, ncol = 2) + geom_hline(yintercept = 0, col = "red") +
    xlab("") + ylab("") + scale_x_continuous(breaks = seq(0, 168, 24)) + theme(legend.position = "none")

2 Clustering sur les coordonnées de l’ACP

Dans cette partie, nous allons utiliser les coordonnées de l’ACP comme matrice de données pour classer les stations.

velibacp <- respca$ind$coord[, 1:7]

2.1 Clustering avec l’algorithme des Kmeans

Question : A l’aide d’une procédure des Kmeans, déterminez une classification des stations.

# A COMPLETER

Correction :

set.seed(12345)
Kmax <- 20
reskmeans <- matrix(0, nrow = nrow(velibacp), ncol = (Kmax - 1))
Iintra <- NULL
Silhou <- NULL
for (k in 2:Kmax) {
    resaux <- kmeans(velibacp, k, nstart = 20, iter.max = 30)
    reskmeans[, (k - 1)] <- resaux$cluster
    Iintra <- c(Iintra, resaux$tot.withinss)
    aux <- silhouette(resaux$cluster, daisy(velibacp))
    Silhou <- c(Silhou, mean(aux[, 3]))
}

rm(resaux, aux)

df <- data.frame(K = 2:Kmax, Iintra = Iintra, Silhou = Silhou)
g1 <- ggplot(df, aes(x = K, y = Iintra)) + geom_line() + geom_point() + xlab("Nombre de classes") +
    ylab("Inertie intraclasse")
g2 <- ggplot(df, aes(x = K, y = Silhou)) + geom_line() + geom_point() + xlab("Nombre de classes") +
    ylab("Critère Silhouette")
grid.arrange(g1, g2, ncol = 2)

resICL <- mclustICL(velibacp, G = 1:20, modelNames = c("EII"))
summary(resICL)

Best ICL values:
            EII,11      EII,13       EII,19
ICL      -41830.01 -41855.1410 -41862.86314
ICL diff      0.00    -25.1267    -32.84888

Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue.

Correction :

On retient la classification en \(K=5\) classes :

Kkmeans <- 5
classifkmeans <- reskmeans[, Kkmeans - 1]
table(classifkmeans)

classifkmeans
  1   2   3   4   5 
406 151 209 241 182 

On visualise la classification sur la carte 2D :

qmplot(longitude, latitude, data = station, colour = as.factor(classifkmeans))

Profil moyen par classe :

plotprofilsmoy(Velibdata, clustering = classifkmeans, plotsem = TRUE)

plotprofilsmoy(Velibdata, clustering = classifkmeans, plotsem = FALSE)

Croisement avec les informations auxilaires sur les stations :

stationcaract(station, clustering = classifkmeans)

Visualisation d’une classe par exemple

I <- which(classifkmeans == 5)
qmplot(longitude, latitude, data = station[I, ])

2.2 Clustering avec de la classification hiérarchique

Question : A l’aide d’une procédure hiérarchique, déterminez une classification des stations.

# A COMPLETER

Correction :

Kmax <- 20
resward <- hclust(dist(velibacp, method = "euclidean"), method = "ward.D2")
df <- data.frame(K = 1:Kmax, height = sort(resward$height, decreasing = T)[1:Kmax])
ggplot(df, aes(x = K, y = height)) + geom_line() + geom_point()

CH <- NULL
Kmax <- 20
for (k in 2:Kmax) {
    CH <- c(CH, index.G1(velibacp, cutree(resward, k)))
}
daux <- data.frame(NbClust = 2:Kmax, CH = CH)
ggplot(daux, aes(x = NbClust, y = CH)) + geom_line() + geom_point()

Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par CAH. Comparez cette classification avec celle obtenue par les Kmeans.

# A COMPLETER

Correction :

Les critères précédents proposent de retenir 3 ou 5 classes :

classifCAH3 <- cutree(resward, 3)
classifCAH5 <- cutree(resward, 5)
table(classifCAH3, classifCAH5)

           classifCAH5
classifCAH3   1   2   3   4   5
          1 163   0   0   0 180
          2   0 237 180   0   0
          3   0   0   0 429   0

On visualise la classification sur la carte 2D :

qmplot(longitude, latitude, data = station, colour = as.factor(classifCAH5))

Profil moyen par classe :

plotprofilsmoy(Velibdata, clustering = classifCAH5, plotsem = TRUE)

plotprofilsmoy(Velibdata, clustering = classifCAH5, plotsem = FALSE)

Croisement avec les informations auxilaires sur les stations :

stationcaract(station, clustering = classifCAH5)

Visualisation d’une classe par exemple

I <- which(classifCAH5 == 2)
qmplot(longitude, latitude, data = station[I, ])

Comparaison avec la classification des Kmeans :

table(classifkmeans, classifCAH5)

             classifCAH5
classifkmeans   1   2   3   4   5
            1   0  22   3 381   0
            2   0 119  26   0   6
            3  12  37   0   0 160
            4 151  33   0  43  14
            5   0  26 151   5   0

adjustedRandIndex(classifkmeans, classifCAH5)

[1] 0.6411121

2.3 Clustering avec des mélanges gaussiens

Question : A l’aide des mélanges gaussiens, déterminez une classification des stations.

# A COMPLETER

Correction :

• Sélection du meilleur modèle de mélange par le critère ICL :

set.seed(12345)
resICLall <- mclustICL(velibacp, G = 2:20)
summary(resICLall)

Best ICL values:
             VEV,7       VVE,10        VVE,7
ICL      -40227.06 -40271.08062 -40281.27660
ICL diff      0.00    -44.02224    -54.21822

resICL <- Mclust(velibacp, G = 7, modelNames = "VEV")

• du meilleur modèle de mélange par le critère BIC

resBICall <- mclustBIC(velibacp, G = 2:20)
summary(resBICall)

Best BIC values:
         VVE,10       VEI,13       VEE,13
BIC      -39953 -39975.77232 -39981.70915
BIC diff      0    -22.77722    -28.71404

resBIC <- Mclust(velibacp, G = 10, modelNames = "VVE")

Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par modèles de mélange. Comparez cette classification avec celles obtenues précédemment.

# A COMPLETER

Correction :

Aux <- data.frame(label = paste("Cl", resICL$classification, sep = ""), proba = apply(resICL$z,
    1, max))
h1 <- ggplot(Aux, aes(x = label, y = proba)) + geom_boxplot()
h2 <- fviz_cluster(resICL, data = velibacp, ellipse = F, geom = "point") + ggtitle("") +
    theme(legend.position = "none")
grid.arrange(h1, h2, ncol = 2)

plotprofilsmoy(Velibdata, clustering = resICL$classification, plotsem = TRUE)

plotprofilsmoy(Velibdata, clustering = resICL$classification, plotsem = FALSE)

qmplot(longitude, latitude, data = station, colour = as.factor(resICL$classification))

table(resICL$classification)

  1   2   3   4   5   6   7 
331 186  83 207 132 128 122 

Comparaison entre la classification avec ICL et celle des Kmeans

table(resICL$classification, classifkmeans)

   classifkmeans
      1   2   3   4   5
  1  42  89  79  89  32
  2  44   0   0 142   0
  3   2   0   0   0  81
  4 190   0   0   8   9
  5   0   0 130   2   0
  6 128   0   0   0   0
  7   0  62   0   0  60

3 Clustering sur un résumé des données

Dans cette partie, on va s’intéresser à un clustering sur un jeu de données “résumé” au sens suivant : on fait la moyenne des loadings selon les tranches horaires suivantes : 0-7h, 8h-20h et 21h-23h pour chaque jour et chaque station.

3.1 Préparation des données

VelibMean <- matrix(0, nrow = nrow(Velibdata), ncol = 3 * 7)
for (jour in 1:7) {
    VelibMean[, 3 * (jour - 1) + 1] <- apply(Velibdata[, 24 * (jour - 1) + c(1:8)],
        1, mean)
    VelibMean[, 3 * (jour - 1) + 2] <- apply(Velibdata[, 24 * (jour - 1) + c(9:21)],
        1, mean)
    VelibMean[, 3 * (jour - 1) + 3] <- apply(Velibdata[, 24 * (jour - 1) + c(22:24)],
        1, mean)
}
colnames(VelibMean) <- paste(rep(unique(day), each = 3), "-", rep(c("0-7h", "8h-20h",
    "21-23h"), 7), sep = "")

Question : Etudiez les corrélations entre variables du jeu de données VelibMean.

# A COMPLETER

Correction :

corrplot(cor(VelibMean[, ((1:7) * 3) - rep(c(2, 1, 0), each = 7)]), method = "ellipse")

Question : Faites une ACP et commentez.

# A COMPLETER

Correction :

resacpMean <- PCA(VelibMean, scale.unit = T, graph = F, ncp = 10)
# Variance expliquée
fviz_eig(resacpMean)

# Visualisation des variables
fviz_pca_var(resacpMean, axes = c(1, 2), geom = c("arrow", "text"))

fviz_pca_var(resacpMean, axes = c(2, 3), geom = c("arrow", "text"))

corrplot(resacpMean$var$cor[, 1:3], method = "ellipse")

# Projection des individus
fviz_pca_ind(resacpMean, axes = c(1, 2), geom = c("point"), habillage = as.factor(station$colline))

fviz_pca_ind(resacpMean, axes = c(1, 3), geom = c("point"))

3.2 Méthodes de classification

3.2.1 Kmeans

Question : Faites une classification des stations à partir du jeu de données VelibMean à l’aide des Kmeans. Etudiez la classification obtenue.

# A COMPLETER

Correction :

set.seed(12345)
Kmax <- 20
reskmeansbis <- matrix(0, nrow = nrow(VelibMean), ncol = (Kmax - 1))
Iintra <- NULL
for (k in 2:Kmax) {
    aux <- kmeans(VelibMean, k)
    reskmeansbis[, (k - 1)] <- aux$cluster
    Iintra <- c(Iintra, aux$tot.withinss)
}

df <- data.frame(K = 2:Kmax, Iintra = Iintra)
ggplot(df, aes(x = K, y = Iintra)) + geom_line() + geom_point()

Etude de la classification obtenue

classifkmeansbis <- reskmeansbis[, 4]
table(classifkmeansbis)

classifkmeansbis
  1   2   3   4   5 
174 254 200 253 308 

plotprofilsmoy(Velibdata, clustering = classifkmeansbis, plotsem = FALSE)

plotprofilsmoy(Velibdata, clustering = classifkmeansbis, plotsem = TRUE)

qmplot(longitude, latitude, data = station, colour = as.factor(classifkmeansbis))

Comparaison avec la classification obtenue sur les coordonnées de l’ACP :

table(classifkmeans, classifkmeansbis)

             classifkmeansbis
classifkmeans   1   2   3   4   5
            1   1   5 121   0 279
            2   9  98   0  44   0
            3   8   0   0 201   0
            4 156   1  76   8   0
            5   0 150   3   0  29

Comparaison avec les variables colline et tourisme

table(classifkmeansbis, station$colline)

                
classifkmeansbis   0   1
               1 172   2
               2 253   1
               3 149  51
               4 252   1
               5 236  72

table(classifkmeansbis, station$tourisme)

                
classifkmeansbis   0   1
               1 174   0
               2 236  18
               3 199   1
               4 248   5
               5 286  22

3.2.2 Mélanges gaussiens

Question : Faites une classification des stations à partir du jeu de données VelibMean à l’aide des mélanges gaussiens. Etudiez la classification obtenue.

# A COMPLETER

Correction :

set.seed(12345)
resICLallbis <- mclustICL(VelibMean, G = 10:20, modelNames = c("VVE", "VEV", "EVV",
    "VVV"))
summary(resICLallbis)

Best ICL values:
           VVE,15     VVE,16    VVE,14
ICL      22181.01 22079.2176 22070.926
ICL diff     0.00  -101.7971  -110.089

resICLbis <- Mclust(VelibMean, G = 16, modelNames = "VVE")

Aux <- data.frame(label = paste("Cl", resICLbis$classification, sep = ""), proba = apply(resICLbis$z,
    1, max))
h1 <- ggplot(Aux, aes(x = label, y = proba)) + geom_boxplot()
h2 <- fviz_cluster(resICLbis, data = velibacp, ellipse = F, geom = "point") + ggtitle("") +
    theme(legend.position = "none")
grid.arrange(h1, h2, ncol = 2)

table(resICLbis$classification)

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
 74 192  92  83  53 125  54  63  52  34  43  54  37 105  54  74 

plotprofilsmoy(Velibdata, clustering = resICLbis$classification, plotsem = FALSE)

plotprofilsmoy(Velibdata, clustering = resICLbis$classification, plotsem = TRUE)

qmplot(longitude, latitude, data = station, colour = as.factor(resICLbis$classification))

stationcaract(station, resICLbis$classification)

4 Clustering pour des séries temporelles

Les données de Vélib sont des données fonctionnelles (on étudie un phénomène qui évolue au cours du temps). Il existe des méthodes de clustering dédiées aux données fonctionnelles, comme par exemple le package funFEM.

Pour utiliser cette méthode, on commence par considérer la décompostion dans la base de Fourier des courbes.

basis <- create.fourier.basis(c(0, 167), nbasis = 25)
fdobj <- smooth.basis(1:168, t(Velibdata), basis)$fd

La méthode funFEM est une procédure de clustering basée sur des mélanges fonctionnels. Comme pour les mélanges gaussiens, il y a différentes formes disponibles selon les contraintes imposés dans le modèle (le détail est hors programme !). Cette méthode est détaillée dans l’article de Bouveyron et al. (2015).

Question : Executez les commandes suivantes et exploitez les résultats. On considère une classification en 10 classes comme dans Bouveyron et al. (2015).

resfunFEM <- funFEM(fdobj, K = 10, model = "AkjBk", init = "kmeans", lambda = 0,
    disp = TRUE)

>> K = 10 
AkjBk   :    bic = -20459.07 
The best model is AkjBk with K = 10 ( bic = -20459.07 )

df <- data.frame(proba = apply(resfunFEM$P, 1, max), classif = as.factor(apply(resfunFEM$P,
    1, which.max)))
ggplot(df, aes(x = classif, y = proba)) + geom_boxplot()

• Etude de la classification obtenue

classifFEM <- apply(resfunFEM$P, 1, which.max)
plotprofilsmoy(Velibdata, clustering = classifFEM, plotsem = FALSE)

plotprofilsmoy(Velibdata, clustering = classifFEM, plotsem = TRUE)

qmplot(longitude, latitude, data = station, colour = as.factor(classifFEM))

stationcaract(station, classifFEM)

• Comparaison avec les classifications précédentes

table(resICL$classification, classifFEM)

   classifFEM
      1   2   3   4   5   6   7   8   9  10
  1  16  59  30  24  24   7   2  48   0 121
  2  13   3 117   0   0   0  37  15   1   0
  3   0   0   0   0   0  70   0  13   0   0
  4   0   0   1   0   0   1  83  88  34   0
  5  59   1   0  70   0   0   0   0   0   2
  6   0   0   0   0   0   0   4   8 116   0
  7   0   0   0   0  65  34   0   2   0  21

adjustedRandIndex(resICL$classification, classifFEM)

[1] 0.3506168

table(resICLbis$classification, classifFEM)

    classifFEM
       1   2   3   4   5   6   7   8   9  10
  1    0   2   7   0   0   0  21  40   0   4
  2   10   2 103   0   0   0  42  24  10   1
  3    0   0   0  13  63   4   0   1   0  11
  4    0   1   0   0   9  18   0  21   0  34
  5    0   0   0   0   5  43   0   5   0   0
  6   14   4   9  24   1   2   0   3   0  68
  7    0   0   1   0   0   0  50   1   2   0
 [ getOption("max.print") est atteint -- 9 lignes omises ]

adjustedRandIndex(resICLbis$classification, classifFEM)

[1] 0.2955738